{
 "metadata": {
  "signature": "sha256:d9abd3b684432593f91b04e433a8ba70537b122f870b89f1135d945b5e41f6c4"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Matplotlib: plotting images with special values\n",
      "======================================================================\n",
      "\n",
      "Image plotting requires data, a colormap, and a normalization. A common\n",
      "desire is to show missing data or other values in a specified color. The\n",
      "following code shows an example of how to do this.\n",
      "\n",
      "The code creates a new Colormap subclass and a norm subclass.\n",
      "\n",
      "The initialization takes a dictionary of value, color pairs. The data is\n",
      "already assumed to be normalized (except for the sentinels which are\n",
      "preserved). The RGB values at the sentinel values are replaced by the\n",
      "specified colors.\n",
      "\n",
      "The class normalizes the data in the standard way except for one\n",
      "sublety. takes an \"ignore\" argument. The ignored values need to be\n",
      "excluded from the normalization so that they do not skew the results.\n",
      "\n",
      "I use a not particularly wonderful algorithm of explicitly sorting the\n",
      "data and using the first non-sentinel values to define the min and max.\n",
      "This can probably be improved, but for my purposes was easy and\n",
      "sufficient. The data is then normalized including the sentinels.\n",
      "Finally, the sentinels are replaced."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib.colors import Colormap, normalize\n",
      "import matplotlib.numerix as nx\n",
      "from types import IntType, FloatType, ListType\n",
      "\n",
      "class SentinelMap(Colormap):\n",
      "        def __init__(self, cmap, sentinels={}):\n",
      "                # boilerplate stuff\n",
      "                self.N = cmap.N\n",
      "                self.name = 'SentinelMap'\n",
      "                self.cmap = cmap\n",
      "                self.sentinels = sentinels\n",
      "                for rgb in sentinels.values():\n",
      "                        if len(rgb)!=3:\n",
      "                                raise ValueError('sentinel color must be RGB')\n",
      "\n",
      "\n",
      "        def __call__(self, scaledImageData, alpha=1):\n",
      "                # assumes the data is already normalized (ignoring sentinels)\n",
      "                # clip to be on the safe side\n",
      "                rgbaValues = self.cmap(nx.clip(scaledImageData, 0.,1.))\n",
      "\n",
      "                #replace sentinel data with sentinel colors\n",
      "                for sentinel,rgb in self.sentinels.items():\n",
      "                        r,g,b = rgb\n",
      "                        rgbaValues[:,:,0] =  nx.where(scaledImageData==sentinel, r, rgbaValues[:,:,0])\n",
      "                        rgbaValues[:,:,1] =  nx.where(scaledImageData==sentinel, g, rgbaValues[:,:,1])\n",
      "                        rgbaValues[:,:,2] =  nx.where(scaledImageData==sentinel, b, rgbaValues[:,:,2])\n",
      "                        rgbaValues[:,:,3] =  nx.where(scaledImageData==sentinel, alpha, rgbaValues[:,:,3])\n",
      "\n",
      "                return rgbaValues\n",
      "\n",
      "class SentinelNorm(normalize):\n",
      "        \"\"\"\n",
      "        Leave the sentinel unchanged\n",
      "        \"\"\"\n",
      "        def __init__(self, ignore=[], vmin=None, vmax=None):\n",
      "                self.vmin=vmin\n",
      "                self.vmax=vmax\n",
      "\n",
      "                if type(ignore) in [IntType, FloatType]:\n",
      "                        self.ignore = [ignore]\n",
      "                else:\n",
      "                        self.ignore = list(ignore)\n",
      "\n",
      "\n",
      "        def __call__(self, value):\n",
      "\n",
      "                vmin = self.vmin\n",
      "                vmax = self.vmax\n",
      "\n",
      "                if type(value) in [IntType, FloatType]:\n",
      "                        vtype = 'scalar'\n",
      "                        val = array([value])\n",
      "                else:\n",
      "                        vtype = 'array'\n",
      "                        val = nx.asarray(value)\n",
      "\n",
      "                # if both vmin is None and vmax is None, we'll automatically\n",
      "                # norm the data to vmin/vmax of the actual data, so the\n",
      "                # clipping step won't be needed.\n",
      "                if vmin is None and vmax is None:\n",
      "                        needs_clipping = False\n",
      "                else:\n",
      "                        needs_clipping = True\n",
      "\n",
      "                if vmin is None or vmax is None:\n",
      "                        rval = nx.ravel(val)\n",
      "                        #do this if sentinels (values to ignore in data)\n",
      "                        if self.ignore:\n",
      "                                sortValues=nx.sort(rval)\n",
      "                                if vmin is None: \n",
      "                                        # find the lowest non-sentinel value\n",
      "                                        for thisVal in sortValues:\n",
      "                                                if thisVal not in self.ignore:\n",
      "                                                        vmin=thisVal #vmin is the lowest non-sentinel value\n",
      "                                                        break\n",
      "                                        else:\n",
      "                                                vmin=0.\n",
      "                                if vmax is None: \n",
      "                                        for thisVal in sortValues[::-1]:\n",
      "                                                if thisVal not in self.ignore:\n",
      "                                                        vmax=thisVal #vmax is the greatest non-sentinel value\n",
      "                                                        break\n",
      "                                        else:\n",
      "                                                vmax=0.\n",
      "                        else:\n",
      "                                if vmin is None: vmin = min(rval)\n",
      "                                if vmax is None: vmax = max(rval)\n",
      "                if vmin > vmax:\n",
      "                        raise ValueError(\"minvalue must be less than or equal to maxvalue\")\n",
      "                elif vmin==vmax:\n",
      "                        return 0.*value\n",
      "                else:\n",
      "                        if needs_clipping:\n",
      "                                val = nx.clip(val,vmin, vmax)\n",
      "                        result = (1.0/(vmax-vmin))*(val-vmin)\n",
      "\n",
      "                # replace sentinels with original (non-normalized) values\n",
      "                for thisIgnore in self.ignore:\n",
      "                        result = nx.where(val==thisIgnore,thisIgnore,result)\n",
      "\n",
      "                if vtype == 'scalar':\n",
      "                        result = result[0]\n",
      "                return result\n",
      "\n",
      "\n",
      "if __name__==\"__main__\":\n",
      "        import pylab\n",
      "        import matplotlib.colors\n",
      "        n=100\n",
      "\n",
      "        # create a random array\n",
      "        X = nx.mlab.rand(n,n)\n",
      "        cmBase = pylab.cm.jet\n",
      "\n",
      "        # plot it array as an image\n",
      "        pylab.figure(1)\n",
      "        pylab.imshow(X, cmap=cmBase, interpolation='nearest')\n",
      "\n",
      "        # define the sentinels\n",
      "        sentinel1 = -10\n",
      "        sentinel2 = 10\n",
      "\n",
      "        # replace some data with sentinels\n",
      "        X[int(.1*n):int(.2*n), int(.5*n):int(.7*n)]  = sentinel1\n",
      "        X[int(.6*n):int(.8*n), int(.2*n):int(.3*n)]  = sentinel2\n",
      "\n",
      "        # define the colormap and norm\n",
      "        rgb1 = (0.,0.,0.)\n",
      "        rgb2 = (1.,0.,0.)\n",
      "        cmap = SentinelMap(cmBase, sentinels={sentinel1:rgb1,sentinel2:rgb2,})\n",
      "        norm = SentinelNorm(ignore=[sentinel1,sentinel2])\n",
      "\n",
      "        # plot with the modified colormap and norm\n",
      "        pylab.figure(2)\n",
      "        pylab.imshow(X, cmap = cmap, norm=norm, interpolation='nearest')\n",
      "\n",
      "        pylab.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If the preceeding code is run from a prompt, two images are generated. The first is\n",
      "a pristine image of random data. The second image is the data modified by setting\n",
      "some blocks to sentinel values and then plotting the sentinels in specific colors.\n",
      "A sample result is shown below.\n"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}